###########################################
# 	BIAS WILL FIND A WAY			#
#	REPLICATION SCRIPT FOR THE APPENDIX	#
#	AUTHOR: MARTIN BISGAARD			#
#	DATE: 06-03-2015				#
###########################################



#	Note: Be sure load the dataset from the main replication 
#	file ("replication-figures2-5.R") and recode the variables before 
#	running this script.

#SIMPLE FUNCTION FOR CLUSTER ROBUST SEs#

cluster.vcov<-function(model, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  
  #calculate degree of freedom adjustment
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- model$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  
  #calculate the uj's
  uj  <- apply(estfun(model),2, function(x) tapply(x, cluster, sum))
    uj  <- na.omit(uj) #if the cluster var is NA at some levels


  #use sandwich to get the var-covar matrix
  vcovCL <- dfc*sandwich(model, meat=crossprod(uj)/N)
  return(vcovCL)
}







###################################################################
#	APPENDIX B: DISTRIBUTION OF RESPONDENTS ACROSS WAVES		#
#	Figure A1									#

pid.dist<-matrix(NA,nc=2,nr=71)
for (i in 1:71) pid.dist[i,]<-table(newdat$party.id2[newdat$date==i])

par(mfrow=c(1,2))
plot(density(pid.dist[,1]),ylim=c(0,.015),xlab="Number of Respondents",main="Labour")
abline(v=mean(pid.dist[,1]),lwd=3)
plot(density(pid.dist[,2]),ylim=c(0,.015),xlab="Number of Respondents",main="Conservative")
abline(v=mean(pid.dist[,2]),lwd=3)


###########################################
# 	APPENDIX C: COVARIATE BALANCE 	#
#	Figure A2					#

varnam<-c("coneu","party.str","mortgage","education","income","gender","age","polatt","unemp","readnwsp")

cobalance.lab<-matrix(NA,nc=length(varnam),nr=71)
cobalance.con<-matrix(NA,nc=length(varnam),nr=71)

for (k in 1:length(varnam)){
for (i in 1:71){
cobalance.lab[i,k]<-mean(newdat[newdat$date==i&newdat$party.id2==0,varnam[k]],na.rm=TRUE)
cobalance.con[i,k]<-mean(newdat[newdat$date==i&newdat$party.id2==1,varnam[k]],na.rm=TRUE)
	}
	}
colnames(cobalance.lab)<-varnam;colnames(cobalance.con)<-varnam

col1<-rgb(0,0,0,alpha=.7)

vartitle<-c("Against British EU Membership","Strength of PID","Has Mortgage","Education","Income","Gender","Year of Birth (0=1900)","Political Attention","Unemployed","Newspaper Exposure")


#pdf('cobalance.pdf',width=8,height=11)

par(mfrow=c(4,3),omd=c(0,1,0,1))
for(k in 1:length(varnam)){
plot(cobalance.lab[,varnam[k]],pch=21,bg=col1,ylim=c(range(newdat[,varnam[k]],na.rm=TRUE)),
	main=vartitle[k],ylab="",xaxt="n",xlab="")
ifelse(k==1,legend('bottomright',legend=c('Labour Identifier','Conservative Identifier'),pch=c(21,1),pt.bg=c(col1,NA),bty='n'),"LOL")

segments(0,max(newdat[,varnam[k]],na.rm=TRUE),0,max(newdat[,varnam[k]],na.rm=TRUE)-abs(sd(newdat[,varnam[k]],na.rm=TRUE)),lwd=3,col="darkgrey")
text(x=5,y=max(newdat[,varnam[k]],na.rm=TRUE),labels="1 sd.",xpd=TRUE,adj=0,srt=270)

axis(1,c(1,10,21,32,44,56,68,79,91,97),labels=c("",2005,2006,2007,2008,2009,2010,2011,2012,""),cex.axis=1)
points(cobalance.con[,varnam[k]],pch=1,col=col1)
}

#dev.off()





###############################################################################
#	APPENDIX D.1: ADDRESSING THE VALIDITY OF THE ATTRIBUTION MEASURE (Q28)	#


#new coding scheme was applied in the CMS for rp measures:
newdat[newdat$date %in% c(66,70:97),]$rpbrown<-recode(newdat[newdat$date %in% c(66,70:97),]$rpbrown,"2=0;3=NA;4=NA;9=NA")
table(newdat$rpbrown);table(newdat$date,newdat$rpbrown)

newdat[newdat$date %in% c(66,70:97),]$rpgov<-recode(newdat[newdat$date %in% c(66,70:97),]$rpgov,"2=0;3=NA;4=NA;9=NA")
table(newdat$rpgov);table(newdat$date,newdat$rpgov)

newdat[newdat$date %in% c(66,70:97),]$rpbush<-recode(newdat[newdat$date %in% c(66,70:97),]$rpbush,"2=0;3=NA;4=NA;9=NA")
table(newdat$rpbush);table(newdat$date,newdat$rpbush)

newdat[newdat$date %in% c(66,70:97),]$rpeu<-recode(newdat[newdat$date %in% c(66,70:97),]$rpeu,"2=0;3=NA;4=NA;9=NA")
table(newdat$rpeu);table(newdat$date,newdat$rpeu)

newdat[newdat$date %in% c(66,70:97),]$rpusgov<-recode(newdat[newdat$date %in% c(66,70:97),]$rpusgov,"2=0;3=NA;4=NA;9=NA")
table(newdat$rpusgov);table(newdat$date,newdat$rpusgov)

newdat[newdat$date %in% c(66,70:97),]$rpukbanks<-recode(newdat[newdat$date %in% c(66,70:97),]$rpukbanks,"2=0;3=NA;4=NA;9=NA")
table(newdat$rpukbanks);table(newdat$date,newdat$rpukbanks)

newdat[newdat$date %in% c(66,70:97),]$rpusbanks<-recode(newdat[newdat$date %in% c(66,70:97),]$rpusbanks,"2=0;3=NA;4=NA;9=NA")
table(newdat$rpusbanks);table(newdat$date,newdat$rpusbanks)

newdat[newdat$date %in% c(66,70:97),]$rppmort<-recode(newdat[newdat$date %in% c(66,70:97),]$rppmort,"2=0;3=NA;4=NA;9=NA")
table(newdat$rppmort);table(newdat$date,newdat$rppmort)

newdat[newdat$date %in% c(66,70:97),]$rpoilco<-recode(newdat[newdat$date %in% c(66,70:97),]$rpoilco,"2=0;3=NA;4=NA;9=NA")
table(newdat$rpoilco);table(newdat$date,newdat$rpoilco)

newdat[newdat$date %in% c(66,70:97),]$rpintlfinance<-recode(newdat[newdat$date %in% c(66,70:97),]$rpintlfinance,"2=0;3=NA;4=NA;9=NA")
table(newdat$rpintlfinance);table(newdat$date,newdat$rpintlfinance)

newdat[newdat$date %in% c(66,70:97),]$rpother<-recode(newdat[newdat$date %in% c(66,70:97),]$rpother,"2=0;3=NA;4=NA;9=NA")
table(newdat$rpother);table(newdat$date,newdat$rpother)

newdat[newdat$date %in% c(66,70:97),]$rpnone<-recode(newdat[newdat$date %in% c(66,70:97),]$rpnone,"2=0;3=NA;4=NA;9=NA")
table(newdat$rpnone);table(newdat$date,newdat$rpnone)



###############################
#	CODE FOR FIGURE A3 	#

dvnames<-c("rpbrown","rpgov","rpbush","rpeu","rpusgov","rpukbanks","rpusbanks","rppmort","rpoilco","rpintlfinance","rpother","rpnone")
length(dvnames);dvnames

boot=1000 #WARNING: Runs slow, use 100 alternatively

ame.out<-as.list(rep(NA,12));names(ame.out)<-dvnames
app.out<-as.list(rep(NA,12));names(ame.out)<-dvnames
tabm<-list(NA)
sumstat<-matrix(NA,nr=4,nc=12)

for (i in 1:12){

model<-substitute(sub~party.id2+party.str+polatt+gender+age+income+mortgage+unemp+education+readnwsp+coneu,list(sub=as.name(dvnames[i])))
m1<-glm(model,data=newdat[newdat$date %in% 1:71,],family=binomial(link="probit"))

tabm[[i]]<-m1

sumstat[,i]<-mtable(m1)$summaries[1:4]

	coefs.out<-matrix(NA,nr=1,nc=2);colnames(coefs.out)<-c("AME","SE")
	coefs1.out<-matrix(NA,nr=1,nc=4);colnames(coefs1.out)<-c("LAB","SE","CONS","SE")

	samp<-mvrnorm(n=boot,mu=coef(m1),Sigma=cluster.vcov(m1,cluster=newdat$date[rownames(newdat)%in%rownames(m1$model)])) #Parametric sim: Resample estimates from a multivariate normal 
					
	#setting data (party id 0/1)
	datx0<-newdat;datx0$party.id2<-0
	datx1<-newdat;datx1$party.id2<-1

	ame.samp1<-rep(NA,boot)
	app.samp<-matrix(c(rep(NA,boot),rep(NA,boot)),nr=boot)

	for (k in 1:boot){
		pred<-NA
		pdat<-NA
	
		m1$coefficients<-samp[k,]
	
		ame.samp1[k]<-mean(predict(m1,newdata=datx1,type="response")-predict(m1,newdata=datx0,type="response"),na.rm=T) #calculate AME for each individual

			pred<-predict(m1,newdata=newdat,type="response")
			pdat<-cbind(newdat[rownames(newdat) %in% names(pred),],pred)
		app.samp[k,]<-c(mean(pdat$pred[pdat$party.id2==0],na.rm=TRUE),mean(pdat$pred[pdat$party.id2==1],na.rm=TRUE))
		
	}
		coefs.out[1,1]<-mean(ame.samp1)
		coefs.out[1,2]<-sd(ame.samp1)

		coefs1.out[1,1]<-mean(app.samp[,1]);coefs1.out[1,2]<-sd(app.samp[,1])
		coefs1.out[1,3]<-mean(app.samp[,2]);coefs1.out[1,4]<-sd(app.samp[,2])


ame.out[[i]]<-coefs.out
app.out[[i]]<-coefs1.out	
	}
rownames(sumstat)<-rownames(mtable(m1)$summaries)[1:4]
ame.mat<-matrix(unlist(ame.out),nr=12,byrow=TRUE);rownames(ame.mat)<-dvnames;colnames(ame.mat)<-c("AME","SE")
app.mat<-matrix(unlist(app.out),nr=12,byrow=TRUE);rownames(app.mat)<-dvnames;colnames(app.mat)<-c("LAB","SE","CONS","SE")


# PLOTTING FIGURE A3#

col1<-rgb(.2,.2,.2,alpha=.8)
col2<-rgb(.4,.4,.4,alpha=.8)

app.mat<-rbind(NA,NA,app.mat)
ame.mat<-rbind(NA,NA,ame.mat)#to preserve top space

jit=.05

label=c("Gordon Brown","British Government","G. W. Bush","European Union","US Government","British Banks","US Banks",
	"People w/ mortgages","Oil companies","International finance","Other","None")

#pdf('rpcrisis.pdf',width=8,height=4)

#PANEL A#
par(mfrow=c(1,1),omd=c(.1,.6,0,1))

plot(app.mat[,1],14:1-jit,pch=21,xaxt="n",yaxt="n",ylab="",xlab="Avg. Predicted Probability",xlim=c(0,1),col="white")

abline(h=1:12,lwd=1.5,lty=1,col="grey")
points(app.mat[,1],14:1-jit,xpd=TRUE,col=col1,bg=col1,pch=21)
points(app.mat[,3],14:1+jit,xpd=TRUE,col=col1)


segments(app.mat[,1],14:1-jit,app.mat[,1]+2*app.mat[,2],14:1-jit,lwd=2);segments(app.mat[,1],14:1-jit,app.mat[,1]-2*app.mat[,2],14:1-jit,lwd=2)
segments(app.mat[,3],14:1+jit,app.mat[,3]-2*app.mat[,4],14:1+jit,lwd=2);segments(app.mat[,3],14:1+jit,app.mat[,3]+2*app.mat[,4],14:1+jit,lwd=2)

abline(h=13)
mtext(3,line=-1,text="Within Group",font=2)
axis(1)
text(app.mat[3,c(1,3)],y=12.75,labels=c("Labour","Conservative"),pos=1,cex=.7,offset=-.05)

par(omd=c(0,1,0,1),new=TRUE)
text(x=-.05, y=12:1, label,xpd=TRUE,pos=2,xpd=TRUE,cex=.8)

#PANEL B#
par(omd=c(.5,1,0,1),new=TRUE)
plot(ame.mat[,1],14:1,pch=21,xaxt="n",yaxt="n",ylab="",xlab="Avg. Marginal Effect",xlim=c(-.1,.4),col="white")

abline(h=1:12,lwd=1.5,lty=1,col="grey")
points(ame.mat[,1],14:1,xpd=TRUE,col=col1,bg=col1,pch=21)

segments(ame.mat[,1],14:1,ame.mat[,1]+1.95*ame.mat[,2],14:1,lwd=2,col=col1)
segments(ame.mat[,1],14:1,ame.mat[,1]-1.95*ame.mat[,2],14:1,lwd=2,col=col1)

abline(h=13)
mtext(3,line=-1,text="Partisan Differences",font=2)
segments(0,0,0,13,lty=2)
axis(1)

#dev.off()


###############################
#	FULL MODEL OUTPUT		#
#	TABLE A2			#
require(stargazer)

covlab<-c("Conservative ID","Strong identifier","Political attention","Gender","Age","Income","Has mortgage","Unemployed","Education","Newspaper exposure","Against EU membership","Constant")
deplab1<-c("Gordon Brown","British Gov.","George Bush","EU","American Gov.","British Banks")
deplab2<-c("American banks","People w/ mortgages","Oil Companies","Int'l Finance","Other Actor","None")
dvnames

stargazer(tabm[1:6],no.space=TRUE,covariate.labels=covlab,dep.var.labels=deplab1,font.size="scriptsize",label="rpcrisis-table",title="The influence of party identification on blaming different actors for the financial crisis. Estimates are probit coefficients.")
stargazer(tabm[7:12],no.space=TRUE,covariate.labels=covlab,dep.var.labels=deplab2,font.size="scriptsize",title="Continued...")

#getting Rsqrts
toLatex(sumstat[,1:6])
toLatex(sumstat[,7:12])







#################################################################
#  APPENDIX D.2: CONVERGENT VALIDITY OF THE ATTRIBUTION MEASURE #
# TABLE A3 & FIGURE A4


fitgov<-glm(rpgov~rpecon2,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71&!is.na(newdat$party.id2),])
fitbrown<-glm(rpbrown~rpecon2,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71&!is.na(newdat$party.id2),])
fiteu<-glm(rpeu~rpecon2,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71&!is.na(newdat$party.id2),])
fitintl<-glm(rpintlfinance~rpecon2,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71&!is.na(newdat$party.id2),])
fitukbanks<-glm(rpukbanks~rpecon2,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71&!is.na(newdat$party.id2),])
fitusbanks<-glm(rpusbanks~rpecon2,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71&!is.na(newdat$party.id2),])
fitusgov<-glm(rpusgov~rpecon2,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71&!is.na(newdat$party.id2),])
fitbush<-glm(rpbush~rpecon2,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71&!is.na(newdat$party.id2),])
fitother<-glm(rpother~rpecon2,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71&!is.na(newdat$party.id2),])
fitmort<-glm(rppmort~rpecon2,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71&!is.na(newdat$party.id2),])
fitnone<-glm(rpnone~rpecon2,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71&!is.na(newdat$party.id2),])
fitoil<-glm(rpoilco~rpecon2,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71&!is.na(newdat$party.id2),])


mtable(fitbrown,fitgov,fiteu,fitbush,fitusgov,fitintl) #first half of TABLE A3
mtable(fitukbanks,fitusbanks,fitmort,fitoil,fitother,fitnone)#second half



#################################
#PREDICTED PROBS AND UNCERTAINTY#


nm<-paste(names(newdat)[grep("^rp",names(newdat))[c(-1,-14:-19)]],"rpecon2",sep="~")
require(R.utils)
nm<-insert(nm,5,nm[2])
nm<-nm[-2]
mlist<-as.list(nm) #set model list for loop


simdat<-data.frame("(Intercept)"=1,"rpecon2"=c(0,1)) #simple predict data

boot=100 #set no. of samples
eff.out<-NA
se.out<-NA
for( i in 1:length(mlist)){
	model<-mlist[[i]]
	fit<-glm(model,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71&!is.na(newdat$party.id2),])
	
	vmat<-vcov(fit,cluster=newdat$date[rownames(newdat)%in%rownames(fit$model)]) #remember to load script
	samp<-mvrnorm(n=boot,mu=coef(fit),Sigma=vmat)

	pred.out<-apply(samp,1,function(x){
	pnorm(as.matrix(simdat)%*%as.matrix(x)) #calculates predicted probs from sampled estimates
	})
	
	ame.out<-apply(pred.out,2,diff) #adding AMEs
	eff.out[i]<-mean(ame.out)
	se.out[i]<-sd(ame.out)
}

eff.out<-rev(eff.out)  #rev for plotting horiz
se.out<-rev(se.out)

###PLOTTING####

col1<-rgb(.2,.2,.2,alpha=.8)
col2<-rgb(.4,.4,.4,alpha=.8)

varlab<-c("Gordon Brown","British government","EU","G. W. Bush","US government","British banks","US banks","People w/ mortgages","Oil companies","International finance","Other","None")

#pdf('convergentvalid.pdf',width=6.75,height=4.5)

par(mfrow=c(1,1),omd=c(.2,1,0,1))
plot(eff.out,1:12,xlim=c(-.2,.3),yaxt="n",ylab="",xlab="")
mtext(1,text="Marginal Effect (Change in Probability)",line=2.25)
mtext(1,text="P(y=1| British Gov. affects economy) - P(y=1| EU or Neither affects economy)",line=3.25,cex=.8)
abline(h=1:12,col="grey",lwd=1.5)
abline(v=0,lty=2)
axis(2,at=1:12,labels=rev(varlab),las=2,cex.axis=1)
segments(eff.out,1:12,eff.out+2*se.out,1:12,col=col1,lwd=1.5)
segments(eff.out,1:12,eff.out-2*se.out,1:12,col=col1,lwd=1.5)
points(eff.out,1:12,pch=21,bg=col1,cex=1)

par(omd=c(0,1,0,1),new=T)
mtext(3,adj=-0.2,text="Who is responsible for",line=1.5,font=3)
mtext(3,adj=-0.12,text="the financial crisis?",line=.5,font=3)
segments(-.3,12.85,-.125,12.85,xpd=T,col="black",lwd=1.5)

#dev.off()




#############################################################
#	APPENDIX D.3: DO DIFFERNET RESPONDENTS ATTRIBUTE 	#
#	DIFFERENT MEANING TO THE RESPONSE CATEGORIES		#
#	TABLE A4

newdat$age.bin<-cut(newdat$age,c(30,40,50,60,70,80,90),include.lowest=T,include.highest=T) #binning age variable
table(newdat$age.bin)
newdat$age.bin<-as.numeric(newdat$age.bin)



fitgov0<-glm(rpgov~rpecon2+party.id2+party.str+natretro+polatt+gender+mortgage+unemp+education+income+readnwsp+age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov<-glm(rpgov~rpecon2*party.id2+rpecon2*party.str+rpecon2*natretro+rpecon2*polatt+rpecon2*gender+rpecon2*mortgage+rpecon2*unemp+rpecon2*education+rpecon2*income+rpecon2*readnwsp+rpecon2*age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov1<-glm(rpgov~rpecon2*party.id2+party.str+natretro+polatt+gender+mortgage+unemp+education+income+readnwsp+age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov2<-glm(rpgov~rpecon2+party.id2+rpecon2*party.str+natretro+polatt+gender+mortgage+unemp+education+income+readnwsp+age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov3<-glm(rpgov~rpecon2+party.id2+party.str+rpecon2*natretro+polatt+gender+mortgage+unemp+education+income+readnwsp+age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov4<-glm(rpgov~rpecon2+party.id2+party.str+natretro+rpecon2*polatt+gender+mortgage+unemp+education+income+readnwsp+age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov5<-glm(rpgov~rpecon2+party.id2+party.str+natretro+polatt+rpecon2*gender+mortgage+unemp+education+income+readnwsp+age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov6<-glm(rpgov~rpecon2+party.id2+party.str+natretro+polatt+gender+rpecon2*mortgage+unemp+education+income+readnwsp+age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov7<-glm(rpgov~rpecon2+party.id2+party.str+natretro+polatt+gender+mortgage+rpecon2*unemp+education+income+readnwsp+age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov8<-glm(rpgov~rpecon2+party.id2+party.str+natretro+polatt+gender+mortgage+unemp+rpecon2*education+income+readnwsp+age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov9<-glm(rpgov~rpecon2+party.id2+party.str+natretro+polatt+gender+mortgage+unemp+education+rpecon2*income+readnwsp+age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov10<-glm(rpgov~rpecon2+party.id2+party.str+natretro+polatt+gender+mortgage+unemp+education+income+rpecon2*readnwsp+age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov11<-glm(rpgov~rpecon2+party.id2+party.str+natretro+polatt+gender+mortgage+unemp+education+income+readnwsp+rpecon2*age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])

mtable(fitgov0,fitgov1,fitgov2,fitgov3,fitgov4,fitgov5,fitgov6,fitgov7,fitgov8,fitgov9,fitgov10,fitgov11)

#F-test: does including all the interaction terms significantly improve the fit of the model?
anova(fitgov0,fitgov,test="F")
lrtest(fitgov0,fitgov)#gives same answer



fitgov0<-glm(rpbrown~rpecon2+party.id2+party.str+natretro+polatt+gender+mortgage+unemp+education+income+readnwsp+age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov<-glm(rpbrown~rpecon2*party.id2+rpecon2*party.str+rpecon2*natretro+rpecon2*polatt+rpecon2*gender+rpecon2*mortgage+rpecon2*unemp+rpecon2*education+rpecon2*income+rpecon2*readnwsp+rpecon2*age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov1<-glm(rpbrown~rpecon2*party.id2+party.str+natretro+polatt+gender+mortgage+unemp+education+income+readnwsp+age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov2<-glm(rpbrown~rpecon2+party.id2+rpecon2*party.str+natretro+polatt+gender+mortgage+unemp+education+income+readnwsp+age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov3<-glm(rpbrown~rpecon2+party.id2+party.str+rpecon2*natretro+polatt+gender+mortgage+unemp+education+income+readnwsp+age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov4<-glm(rpbrown~rpecon2+party.id2+party.str+natretro+rpecon2*polatt+gender+mortgage+unemp+education+income+readnwsp+age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov5<-glm(rpbrown~rpecon2+party.id2+party.str+natretro+polatt+rpecon2*gender+mortgage+unemp+education+income+readnwsp+age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov6<-glm(rpbrown~rpecon2+party.id2+party.str+natretro+polatt+gender+rpecon2*mortgage+unemp+education+income+readnwsp+age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov7<-glm(rpbrown~rpecon2+party.id2+party.str+natretro+polatt+gender+mortgage+rpecon2*unemp+education+income+readnwsp+age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov8<-glm(rpbrown~rpecon2+party.id2+party.str+natretro+polatt+gender+mortgage+unemp+rpecon2*education+income+readnwsp+age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov9<-glm(rpbrown~rpecon2+party.id2+party.str+natretro+polatt+gender+mortgage+unemp+education+rpecon2*income+readnwsp+age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov10<-glm(rpbrown~rpecon2+party.id2+party.str+natretro+polatt+gender+mortgage+unemp+education+income+rpecon2*readnwsp+age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])
fitgov11<-glm(rpbrown~rpecon2+party.id2+party.str+natretro+polatt+gender+mortgage+unemp+education+income+readnwsp+rpecon2*age.bin,family=binomial(link="probit"),data=newdat[newdat$date%in%1:71,])

#F-test
anova(fitgov0,fitgov,test="F")
lrtest(fitgov0,fitgov)

mtable(fitgov0,fitgov1,fitgov2,fitgov3,fitgov4,fitgov5,fitgov6,fitgov7,fitgov8,fitgov9,fitgov10,fitgov11)







#################################################
#	E: ADDRESSING POTENTIAL SELECTION BIAS	#
#	FIGURE A5


propout<-matrix(NA,nr=71,nc=2)
labcount<-rep(NA,71)
concount<-rep(NA,71)
for (i in 1:71){
	
	labcount[i]<-length(na.exclude(newdat$party.id2[newdat$party.id2==0 & newdat$date==i]))
	concount[i]<-length(na.exclude(newdat$party.id2[newdat$party.id2==1 & newdat$date==i]))
	samptot<-dim(newdat[newdat$date==i,])[1]	
	propout[i,1]<-labcount[i]/samptot
	propout[i,2]<-concount[i]/samptot
}


#pdf('idgraph.pdf',width=8,height=6)

col1=rgb(.3,.3,.3,alpha=.8)
par(mfrow=c(1,1),omd=c(0,1,0,1))

plot(propout[,1],type="l",lwd=2,col="white",ylim=c(0,.6),xaxt="n",xlab="",ylab="Pct. of total sample")
points(propout[,1],pch=21,bg=col1,col=col1,cex=1)
points(propout[,2],pch=1,col=col1,cex=1,lwd=1)
lines(propout[,2],lwd=1,col="black")
lines(propout[,1],lwd=1,col="black")
axis(1,c(1,10,21,32,44,56,68,79,91,97),labels=c("",2005,2006,2007,2008,2009,2010,2011,2012,""),cex.axis=1)
abline(v=40,lty=2)
text(y=.5,x=42,labels="Northern Rock",cex=1,xpd=TRUE,font=2,srt=270)
legend("topleft",bty="n",pch=c(21,1),pt.bg=c(col1,NA),col=col1,legend=c("Proportion of Labour Id.","Proportion of Conservative Id."),cex=1,lwd=1)

#dev.off()




###############################################################################
#	E.1: ADDRESSING THE ISSUE OF REVERSE CAUSATION USING THE NINE-WAVE	#
#	BES PANEL DATA SET (2005-2010)							#

#First bit of code is for cleaning and reshaping the panel data

load("paneldat.wide.RData")

dim(paneldat.wide)#should be 7793 by 46
names(paneldat.wide)

##### CODING TIMEINVAR PID #####

paneldat.wide$pid.invar<-paneldat.wide$pid.1
paneldat.wide$econ.invar<-paneldat.wide$econ.1
table(paneldat.wide$pid.invar)


#RESHAPING INTO LONG
paneldat<-reshape(paneldat.wide,direction="long",idvar="id",sep=".",varying=c(1:27))
dim(paneldat)
paneldat[1:30,]
names(paneldat)

paneldat<-paneldat[order(paneldat$id),]


############
# RECODING #

paneldat$econ<-recode(paneldat$econ,"6=NA");paneldat$econ<-(paneldat$econ-1)/4
table(paneldat$econ)

paneldat$pid<-recode(paneldat$pid,'1=0;2=1;else=NA')
paneldat$pid.invar<-recode(paneldat$pid.invar,'1=0;2=1;else=NA')


############################################
#	SIMPLE TRANSITION & MISSINGNESS	 #
#	FIGURE A6					#


#first cleaning up wide paneldata
paneldat.wide$econ.1<-recode(paneldat.wide$econ.1,'6=NA')
paneldat.wide$econ.pro<-recode(paneldat.wide$econ.pro,'6=NA')
paneldat.wide$pid.strength<-recode(paneldat.wide$pid.strength,'4=NA')
paneldat.wide$resp.1<-recode(paneldat.wide$resp.1,'5=NA')
paneldat.wide$age<-recode(paneldat.wide$age,'1:5=NA')
paneldat.wide$pid.invar<-recode(paneldat.wide$pid.invar,'1=0;2=1;else=NA')

pnam<-paste("pid.",1:9,sep="")
for (i in 1:9) paneldat.wide[,pnam[i]]<-recode(paneldat.wide[,pnam[i]],'1=0;2=1;else=NA')
table(paneldat.wide$pid.2)

paneldat.wide$missing<-!complete.cases(paneldat.wide[,c(paste("pid.",1:9,sep=""),paste("econ.",1:9,sep=""))]) #missingness and attrition is treated as the same here

#only looking at LAB and CON respondents who participated in w1
paneldat.wide<-paneldat.wide[!is.na(paneldat.wide$pid.1),]
dim(paneldat.wide)


#PLOTTING FIGURE A6#

xval<-c(3,4,5,18,42,54,62,63,64)

missing.time<-(1-missing.time)
trans<-(1-trans)


#pdf('paneltrans.pdf',width=8,height=4)
par(mfrow=c(1,2))
plot(xval,c(NA,missing.time[,1]),type="l",ylim=c(0,1),xaxt="n",las=2,ylab="Probability",xlab="")
axis(1,c(12,36,60),labels=c(2006,2008,2010))
#axis(1,1,labels="(Ref)",cex.axis=.8,las=1,line=.8,tick=FALSE)
points(xval,c(NA,missing.time[,1]),pch=21,bg="black")
lines(xval,c(NA,missing.time[,2]),lty=2);points(xval,c(NA,missing.time[,2]))
mtext(3,text="A",line=2,cex=1.5)
mtext(3,text=expression("Pr(Missingness) | party identification"["t=1"]),line=.8,cex=1)

plot(xval,c(NA,trans[,1]),type="l",ylim=c(0,1),xaxt="n",,xlab="",ylab="",las=2)
axis(1,c(12,36,60),labels=c(2006,2008,2010))
points(xval,c(NA,trans[,1]),pch=21,bg="black")
lines(xval,c(NA,trans[,2]),lty=2);points(xval,c(NA,trans[,2]))
mtext(3,text="B",line=2,cex=1.5)
mtext(3,text=expression("Pr(Transition) | party identification"["t=1"]),line=.8,cex=1)
legend("topleft",bty="n",pch=c(21,1),pt.bg=c("black",NA),lty=c(1,2),legend=c("Labour Identifier","Conservative identifier"),cex=1)

#dev.off()





############
# TABLE A6 #

require(plm)

names(paneldat)
paneldat<-plm.data(paneldat,indexes=c("id","time"))
pfit.lab<-plm(econ~factor(time)*pid.invar,model="random",data=paneldat)
summary(pfit.lab)

#hausmantest
pfit1<-plm(econ~factor(time)*pid.invar,model="within",data=paneldat)
phtest(pfit1,pfit.lab)#RE model is consistent


vcov.adj<-vcovHC(pfit.lab,cluster="group",type="HC1")

varnam<-c("Constant",paste("Wave",2:9,sep=" "),"Conservative ID",paste("ConID x Wave",2:9,sep=" "))
title<-"Estimates from random intercept model explaining retrospective  perceptions of the national economy conditional on party identification and time of  interview"

stargazer(pfit.lab,se=sqrt(diag(vcov.adj)),p.auto=FALSE,single.row=TRUE,t.auto=FALSE,
	title=title,label="paneltable",dep.var.caption="Economic perceptions",intercept.bottom=FALSE,
	covariate.labels=varnam,
	notes="$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01. Standard errors are clustered with respect to time.")
	



##############
# FIGURE A7 ##

paneldat$pid.invar2<-recode(paneldat$pid.invar,'0=1;1=0') #reverse coding for Conservatives as ref group
pfit.con<-plm(econ~factor(time)*pid.invar2,model="random",data=paneldat)

(pfit.adjL<-coeftest(pfit.lab,vcov=function(x) vcovHC(x,cluster="group",type="HC1")))
(pfit.adjC<-coeftest(pfit.con,vcov=function(x) vcovHC(x,cluster="group",type="HC1")))


coefs<-matrix(c(pfit.adjL[c("(Intercept)",paste("factor(time)",2:9,sep="")),1],
	pfit.adjC[c("(Intercept)",paste("factor(time)",2:9,sep="")),1]),nc=2)
coefs<-rbind(coefs[1,],t(coefs[1,]+t(coefs[2:9,])))


cov1<-(vcovHC(pfit.lab,cluster="group",type="HC1"))
cov2<-(vcovHC(pfit.con,cluster="group",type="HC1"))

se<-matrix(NA,nc=2,nr=8)
for (k in 1:8) se[k,]<-c(sqrt(cov1[k,k]+cov1[k+1,k+1]+2*cov1[k,k+1]),sqrt(cov2[k,k]+cov2[k+1,k+1]+2*cov2[k,k+1]))
se<-rbind(c(sqrt(cov1[1,1]),sqrt(cov1[2,2])),se)

gap<-pfit.lab$coefficients[10:18]
gap<-gap[1]+gap[2:9];gap<-c(pfit.lab$coefficients[10],gap)

gap.se<-NA
for(k in 1:8) gap.se[k]<-sqrt(cov1[k+9,k+9]+cov1[k+10,k+10]+2*cov1[k+9,k+10])
gap.se<-c(sqrt(cov1[10,10]),gap.se)


######## PLOTTING RESULTS########


#pdf('panelmain.pdf',width=8,height=4)

par(mfrow=c(1,2))
plot(xval,coefs[,1],ylim=c(0,1),pch=21,bg="black",xaxt="n",ylab="Retrospective Economic Perceptions",xaxt="n",xlab="")

lines(1:64,SMA(abs(coefs.out[9:72,1]),n=2),lwd=1) #Adding RCS curve from CMS data. NOTE: Run main script for Figure 2
lines(1:64,SMA(abs(coefs.out[9:72,3]),n=2),lwd=1) 

points(xval,coefs[,2],ylim=c(0,1))
segments(xval,coefs[,1],xval,coefs[,1]+2*se[,1]);segments(xval,coefs[,1],xval,coefs[,1]-2*se[,1])
segments(xval,coefs[,2],xval,coefs[,2]+2*se[,2]);segments(xval,coefs[,2],xval,coefs[,2]-2*se[,2])
axis(1,c(12,36,60),labels=c(2006,2008,2010))
mtext(3,text="A",line=2,cex=1.5)
mtext(3,text="Within Group Changes",line=1,cex=1)
legend("topright",bty="n",pch=c(21,1,NA),pt.bg=c("black",NA),lty=c(NA,NA,1),legend=c("Labour identifier","Conservative identifier","Est. CMS curve"),cex=1)


plot(xval,-gap,ylim=c(0,.5),pch=21,bg="grey",xaxt="n",xlab="",ylab="Partisan Difference")
lines(1:64,SMA(abs(coefs.out[9:72,5]),n=2),lwd=1) #Adding RCS curve from CMS data. NOTE: Rune main script before this
segments(xval,-gap,xval,-gap+2*gap.se);segments(xval,-gap,xval,-gap-2*gap.se)
points(xval,-gap,pch=21,bg="grey")
abline(h=0,lty=2)
axis(1,c(12,36,60),labels=c(2006,2008,2010))
mtext(3,text="B",line=2,cex=1.5)
mtext(3,text="Between Group Changes",line=1,cex=1)

#dev.off()






#######################################################################
# APPENDIX G: REPLICATION OF FIGURE 5 WITHOUT COLLAPSING SURVEY WAVES #'


##WEAK ID###
coefs.out<-matrix(NA,nc=6,nr=97)
boot<-100 #set number of reps for AME, WARNING: boot=1000 runs slow
set.seed(123)

lmfit<-glm(rpecon2~party.id2*date+polatt+coneu+gender+age+income+mortgage+unemp+education,
	data=newdat[newdat$party.str==0,],
	family=binomial(link="probit"))
mtable(lmfit)

coefsrp<-matrix(NA,nc=6,nr=71) #matrix to hold overall output for plotting

varnam<-c("coneu","party.str","mortgage","education","income","gender","age","polatt","unemp","party.id2","date")

################
# APP and AME  #

	x<-lmfit 	#Set model
	vmat<-vcov(x) #var-covar matrix
	samp<-mvrnorm(n=boot,mu=coef(x),Sigma=vmat) #simulate estimation uncertainty

	#set up date frame with observered values
	datx1<-newdat[complete.cases(newdat[,varnam]),varnam];datx1$party.id2<-1 #data where all x's are set to 1
	datx0<-newdat[complete.cases(newdat[,varnam]),varnam];datx0$party.id2<-0 #set to 0

	ame.samp<-rep(NA,boot)
	app1<-rep(NA,boot)
	app0<-rep(NA,boot)
	ame<-rep(NA,boot)

	for (i in 1:71){
	for (k in 1:boot){
		x$coefficients<-samp[k,]
		app0[k]<-mean(predict(x,newdat=datx0[datx0$date==i,],type="response"))
		app1[k]<-mean(predict(x,newdat=datx1[datx1$date==i,],type="response"))
		ame[k]<-app0[k]-app1[k]
	}
		coefsrp[i,]<-c(mean(app0),sd(app0),mean(app1),sd(app1),mean(ame),sd(ame)) #storing time varying APPs and AMEs
	}

coefsrp
coefsrp[,5]<- -coefsrp[,5] #reverse gap



#plotting first panel

col1<-rgb(.2,.2,.2,alpha=.8)
col2<-rgb(.4,.4,.4,alpha=.8)

#pdf('rpecon-weak-nocollapse.pdf',width=5.75,height=6.25)

par(mfrow=c(1,1))
par(omd=c(0,1,.375,1))

plot(coefsrp[,1],ylim=c(0,1.15),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="Avg. Predicted Probability",
	col=col1)
mtext(3,line=1.25,text='Weak Identifiers',cex=1.5,font=2)

segments(1:97,coefsrp[,1],1:97,coefsrp[,1]+1.97*coefsrp[,2],col=col1);segments(1:97,coefsrp[,1],1:97,coefsrp[,1]-1.97*coefsrp[,2],col=col1)
segments(1:97,coefsrp[,3],1:97,coefsrp[,3]+1.97*coefsrp[,4],col=col1);segments(1:97,coefsrp[,3],1:97,coefsrp[,3]-1.97*coefsrp[,4],col=col1)

points(coefsrp[,1],pch=21,bg=col1,col=col1)
points(coefsrp[,3],pch=21,col=col1)

lines(SMA(coefsrp[,1],n=2))
lines(SMA(coefsrp[,3],n=2),lty=1)


abline(h=1.05)
mtext(3,text="Attribution of Responsibility",line=-1,font=1)

legend("bottomleft",pch=c(21,1),pt.bg=c(col1,NA),col=col1,legend=c("Labour identifier","Conservative Identifier"),cex=.9,bty="n")

par(omd=c(0,1,.375,.95),new=TRUE)
abline(v=40,lty=2)
par(omd=c(0,1,.375,1),new=TRUE)

text(35,.975,labels="Northern Rock",font=1,cex=1,pos=2)
arrows(33,.975,40,.95,length=.075)

par(omd=c(0,1,0,.625),new=TRUE)



plot(coefsrp[,5],ylim=c(-.5,.75),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="Avg. Marginal Effect",yaxt="n",
	col=col2)
axis(2,at=seq(-.4,.4,.2),las=2)

abline(h=.65)
mtext(3,text="Partisan Differences",line=-1,font=1)

segments(1:97,coefsrp[,5],1:97,coefsrp[,5]+1.97*coefsrp[,6],col=col2);segments(1:97,coefsrp[,5],1:97,coefsrp[,5]-1.97*coefsrp[,6],col=col2)
lines(SMA(coefsrp[,5],n=2))
points(coefsrp[,5],pch=21,bg=col2,col=col2)
abline(h=0,lty=2)

par(omd=c(0,1,0,.575),new=TRUE)
abline(v=40,lty=2)
par(omd=c(0,1,0,.625),new=TRUE)

axis(1,c(10,21,32,44,56,68,79,91),labels=c(2005,"",2007,"",2009,"",2011,""),cex.axis=1)

#dev.off()




###MIDDLING ID####

lmfit<-glm(rpecon2~party.id2*date+polatt+coneu+gender+age+income+mortgage+unemp+education,
	data=newdat[newdat$party.str==.5,],
	family=binomial(link="probit"))

mtable(lmfit)

coefsrp<-matrix(NA,nc=6,nr=71) #matrix to hold overall output for plotting

varnam<-c("coneu","party.str","mortgage","education","income","gender","age","polatt","unemp","party.id2","date")

################
# APP and AME  #

	x<-lmfit 	#Set model
	vmat<-vcov(x) #var-covar matrix
	samp<-mvrnorm(n=boot,mu=coef(x),Sigma=vmat) #simulate estimation uncertainty

	#set up date frame with observered values
	datx1<-newdat[complete.cases(newdat[,varnam]),varnam];datx1$party.id2<-1 #data where all x's are set to 1
	datx0<-newdat[complete.cases(newdat[,varnam]),varnam];datx0$party.id2<-0 #set to 0

	ame.samp<-rep(NA,boot)
	app1<-rep(NA,boot)
	app0<-rep(NA,boot)
	ame<-rep(NA,boot)

	for (i in 1:71){
	for (k in 1:boot){
		x$coefficients<-samp[k,]
		app0[k]<-mean(predict(x,newdat=datx0[datx0$date==i,],type="response"))
		app1[k]<-mean(predict(x,newdat=datx1[datx1$date==i,],type="response"))
		ame[k]<-app0[k]-app1[k]
	}
		coefsrp[i,]<-c(mean(app0),sd(app0),mean(app1),sd(app1),mean(ame),sd(ame)) #storing time varying APPs and AMEs
	}
coefsrp
coefsrp[,5]<- -coefsrp[,5]



#PLOT middle panel#

col1<-rgb(.2,.2,.2,alpha=.8)
col2<-rgb(.4,.4,.4,alpha=.8)

#pdf('rpecon-mid-nocollapse.pdf',width=5.75,height=6.25)

par(mfrow=c(1,1))
par(omd=c(0,1,.375,1))

plot(coefsrp[,1],ylim=c(0,1.15),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="Avg. Predicted Probability",yaxt="n",
	col=col1)
mtext(3,line=1.25,text='Middling Identifiers',cex=1.5,font=2)
axis(2,labels=NA)

segments(1:97,coefsrp[,1],1:97,coefsrp[,1]+1.97*coefsrp[,2],col=col1);segments(1:97,coefsrp[,1],1:97,coefsrp[,1]-1.97*coefsrp[,2],col=col1)
segments(1:97,coefsrp[,3],1:97,coefsrp[,3]+1.97*coefsrp[,4],col=col1);segments(1:97,coefsrp[,3],1:97,coefsrp[,3]-1.97*coefsrp[,4],col=col1)

points(coefsrp[,1],pch=21,bg=col1,col=col1)
points(coefsrp[,3],pch=21,col=col1)

lines(SMA(coefsrp[,1],n=2))
lines(SMA(coefsrp[,3],n=2),lty=1)


abline(h=1.05)
mtext(3,text="Attribution of Responsibility",line=-1,font=1)

legend("bottomleft",pch=c(21,1),pt.bg=c(col1,NA),col=col1,legend=c("Labour identifier","Conservative Identifier"),cex=.9,bty="n")

par(omd=c(0,1,.375,.95),new=TRUE)
abline(v=40,lty=2)
par(omd=c(0,1,.375,1),new=TRUE)

#text(35,.975,labels="Northern Rock",font=1,cex=1,pos=2)
#arrows(33,.975,40,.95,length=.075)

par(omd=c(0,1,0,.625),new=TRUE)



plot(coefsrp[,5],ylim=c(-.5,.75),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="Avg. Marginal Effect",yaxt="n",
	col=col2)
axis(2,at=seq(-.4,.4,.2),labels=NA,las=2)

abline(h=.65)
mtext(3,text="Partisan Differences",line=-1,font=1)

segments(1:97,coefsrp[,5],1:97,coefsrp[,5]+1.97*coefsrp[,6],col=col2);segments(1:97,coefsrp[,5],1:97,coefsrp[,5]-1.97*coefsrp[,6],col=col2)
lines(SMA(coefsrp[,5],n=2))
points(coefsrp[,5],pch=21,bg=col2,col=col2)
abline(h=0,lty=2)

par(omd=c(0,1,0,.575),new=TRUE)
abline(v=40,lty=2)
par(omd=c(0,1,0,.625),new=TRUE)

axis(1,c(10,21,32,44,56,68,79,91),labels=c(2005,"",2007,"",2009,"",2011,""),cex.axis=1)

#dev.off()





##STRONG ID###
lmfit<-glm(rpecon2~party.id2*date+polatt+coneu+gender+age+income+mortgage+unemp+education,
	data=newdat[newdat$party.str==1,],
	family=binomial(link="probit"))

mtable(lmfit)

coefsrp<-matrix(NA,nc=6,nr=71)

varnam<-c("coneu","party.str","mortgage","education","income","gender","age","polatt","unemp","party.id2","date")

################
# APP and AME  #

	x<-lmfit 	#Set model
	vmat<-vcov(x) #var-covar matrix
	samp<-mvrnorm(n=boot,mu=coef(x),Sigma=vmat) #simulate estimation uncertainty

	#set up date frame with observered values
	datx1<-newdat[complete.cases(newdat[,varnam]),varnam];datx1$party.id2<-1 #data where all x's are set to 1
	datx0<-newdat[complete.cases(newdat[,varnam]),varnam];datx0$party.id2<-0 #set to 0

	ame.samp<-rep(NA,boot)
	app1<-rep(NA,boot)
	app0<-rep(NA,boot)
	ame<-rep(NA,boot)

	for (i in 1:71){
	for (k in 1:boot){
		x$coefficients<-samp[k,]
		app0[k]<-mean(predict(x,newdat=datx0[datx0$date==i,],type="response"))
		app1[k]<-mean(predict(x,newdat=datx1[datx1$date==i,],type="response"))
		ame[k]<-app0[k]-app1[k]
	}
		coefsrp[i,]<-c(mean(app0),sd(app0),mean(app1),sd(app1),mean(ame),sd(ame)) #storing time varying APPs and AMEs
	}
coefsrp
coefsrp[,5]<- -coefsrp[,5]



#plotting last panel#
col1<-rgb(.2,.2,.2,alpha=.8)
col2<-rgb(.4,.4,.4,alpha=.8)

#pdf('rpecon-strong-nocollapse.pdf',width=5.75,height=6.25)

par(mfrow=c(1,1))
par(omd=c(0,1,.375,1))
plot(coefsrp[,1],ylim=c(0,1.15),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="",yaxt='n',
	col=col1)
mtext(3,line=1.25,text='Strong Identifiers',cex=1.5,font=2)
axis(2,seq(0,1,.2),labels=NA)
segments(1:97,coefsrp[,1],1:97,coefsrp[,1]+1.97*coefsrp[,2],col=col1);segments(1:97,coefsrp[,1],1:97,coefsrp[,1]-1.97*coefsrp[,2],col=col1)
segments(1:97,coefsrp[,3],1:97,coefsrp[,3]+1.97*coefsrp[,4],col=col1);segments(1:97,coefsrp[,3],1:97,coefsrp[,3]-1.97*coefsrp[,4],col=col1)
points(coefsrp[,1],pch=21,bg=col1,col=col1)
points(coefsrp[,3],pch=21,col=col1)
lines(SMA(coefsrp[,1],n=2))
lines(SMA(coefsrp[,3],n=2),lty=1)
abline(h=1.05)
mtext(3,text="Attribution of Responsibility",line=-1,font=1)
par(omd=c(0,1,.375,.95),new=TRUE)
abline(v=40,lty=2)
par(omd=c(0,1,.375,1),new=TRUE)
par(omd=c(0,1,0,.625),new=TRUE)
plot(coefsrp[,5],ylim=c(-.4,.75),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="",yaxt="n",
	col=col2)
axis(2,at=seq(-.4,.4,.2),labels=NA)
abline(h=.65)
mtext(3,text="Partisan Differences",line=-1,font=1)
segments(1:97,coefsrp[,5],1:97,coefsrp[,5]+1.97*coefsrp[,6],col=col2);segments(1:97,coefsrp[,5],1:97,coefsrp[,5]-1.97*coefsrp[,6],col=col2)
lines(SMA(coefsrp[,5],n=2))
points(coefsrp[,5],pch=21,bg=col2,col=col2)
abline(h=0,lty=2)
par(omd=c(0,1,0,.575),new=TRUE)
abline(v=40,lty=2)
par(omd=c(0,1,0,.625),new=TRUE)
axis(1,c(10,21,32,44,56,68,79,91),labels=c(2005,"",2007,"",2009,"",2011,""),cex.axis=1)

#dev.off()





#################################################
#	APPENDIX H: USING AN ORDERED PROBIT MODEL	#
#	FIGURE A8						#


coefs.outL<-matrix(NA,nc=5,nr=71)
coefs.outC<-matrix(NA,nc=5,nr=71)
simdat<-data.frame("party.id2"=c(0,1))


for (i in 1:71){

plfit<-polr(as.factor(natretro)~party.id2,
	data=newdat[newdat$date==i,],method="probit")
	predval<-predict(plfit,type="probs",newdata=simdat)
	ifelse(length(predval[1,])<5,coefs.outL[i,]<-c(predval[1,],0),coefs.outL[i,]<-predval[1,])
	ifelse(length(predval[2,])<5,coefs.outC[i,]<-c(predval[2,],0),coefs.outC[i,]<-predval[2,])
}


#plotting#	
#pdf('ordered-lab.pdf',height=4,width=6.5)

par(mfrow=c(1,1),omd=c(0,.75,0,1))
plot(coefs.outL[,1],type="line",ylim=c(0,1),xlab="",main="Labour",ylab="Avg. Cumulative Probability",xaxt="n",bty="n")

X.Vec <- c(1:71, rev(1:71),1)
Y.Vec <- c(coefs.outL[,1], rep(0,71),0)
polygon(X.Vec,Y.Vec,col="grey30")

Y.Vec <- c((coefs.outL[,2]+coefs.outL[,1]),rev(coefs.outL[,1]),coefs.outL[1,1])
polygon(X.Vec,Y.Vec,col="grey40")

upper<-coefs.outL[,1]+coefs.outL[,2]
Y.Vec <- c((coefs.outL[,3]+upper),rev(upper),upper[1])
polygon(X.Vec,Y.Vec,col="grey50")

upper<-coefs.outL[,1]+coefs.outL[,2]+coefs.outL[,3]
Y.Vec <- c((coefs.outL[,4]+upper),rev(upper),upper[1])
polygon(X.Vec,Y.Vec,col="grey60")

upper<-coefs.outL[,1]+coefs.outL[,2]+coefs.outL[,3]+coefs.outL[,4]
Y.Vec <- c((coefs.outL[,5]+upper),rev(upper),upper[1])
polygon(X.Vec,Y.Vec,col="grey70")
axis(1,c(1,10,21,32,44,56,68,79,91,97),labels=c("",2005,2006,2007,2008,2009,2010,2011,2012,""),cex.axis=0.8)
abline(v=71,lty=1,lwd=3,col="white")
abline(v=40,lty=1,lwd=3,col="white")
text(y=1.1,x=40,labels="Northern Rock",cex=0.8,xpd=TRUE,font=2)

par(omd=c(.61,1,0,1),new=TRUE)
plot.new()
legend(x=0,y=.9,xpd=TRUE,bty="n",pt.bg=rev(c("grey30","grey40","grey50","grey60","grey70"))
	,pch=22,title=expression(paste("Pr(Y=y"[i], ")")),pt.cex=3.1,
	legend=rev(c(" Got a lot worse"," Got worse",
	" Stayed around the same"," Got better"," Got a lot better")))

#dev.off()

	

#pdf('ordered-con.pdf',height=4,width=6.5)
#CONS

par(mfrow=c(1,1),omd=c(0,.75,0,1))
plot(coefs.outC[,1],type="line",ylim=c(0,1),xlab="",main="Conservative",ylab="Avg. Cumulative Probability",xaxt="n",bty="n")

X.Vec <- c(1:71, rev(1:71),1)
Y.Vec <- c(coefs.outC[,1], rep(0,71),0)
polygon(X.Vec,Y.Vec,col="grey30")

Y.Vec <- c((coefs.outC[,2]+coefs.outC[,1]),rev(coefs.outC[,1]),coefs.outC[1,1])
polygon(X.Vec,Y.Vec,col="grey40")

upper<-coefs.outC[,1]+coefs.outC[,2]
Y.Vec <- c((coefs.outC[,3]+upper),rev(upper),upper[1])
polygon(X.Vec,Y.Vec,col="grey50")

upper<-coefs.outC[,1]+coefs.outC[,2]+coefs.outC[,3]
Y.Vec <- c((coefs.outC[,4]+upper),rev(upper),upper[1])
polygon(X.Vec,Y.Vec,col="grey60")

upper<-coefs.outC[,1]+coefs.outC[,2]+coefs.outC[,3]+coefs.outC[,4]
Y.Vec <- c((coefs.outC[,5]+upper),rev(upper),upper[1])
polygon(X.Vec,Y.Vec,col="grey70")

axis(1,c(1,10,21,32,44,56,68,79,91,97),labels=c("",2005,2006,2007,2008,2009,2010,2011,2012,""),cex.axis=0.8)
abline(v=71,lty=1,lwd=3,col="white")
abline(v=40,lty=1,lwd=3,col="white")
text(y=1.1,x=40,labels="Northern Rock",cex=0.8,xpd=TRUE,font=2)
#text(y=1.1,x=71,labels="Election",cex=0.8,xpd=TRUE,font=2)

par(omd=c(.61,1,0,1),new=TRUE)
plot.new()
legend(x=0,y=.9,xpd=TRUE,bty="n",pt.bg=rev(c("grey30","grey40","grey50","grey60","grey70"))
	,pch=22,title=expression(paste("Pr(Y=y"[i], ")")),pt.cex=3.1,
	legend=rev(c(" Got a lot worse"," Got worse",
	" Stayed around the same"," Got better"," Got a lot better")))

#dev.off()




###############################
#	DIRECT COMPARISONS	#




#PR(GOT WORSE)

fit1<-coefs.outL[42:71,2]+coefs.outL[42:71,1]
fit2<-coefs.outC[42:71,2]+coefs.outC[42:71,1]
dif<-fit1-fit2



par(mfrow=c(2,2))

plot(0:29,fit1,ylim=c(0,1),type="l",
	main=expression("Pr("<= "Got worse | PartyID)"),bty="n",col="white",
	cex.lab=1,ylab="Average Predicted Probability",xlab="Month post outbreak of crisis")

lines(0:29,fit1,col="indianred",lwd=2)
points(0:29,fit1,bg="indianred",pch=21)
lines(0:29,fit2,col="skyblue3",lwd=2)
points(0:29,fit2,bg="skyblue3",pch=21)


legend("bottomleft",legend=c("Labour","Conservative"),pch=21,
	pt.cex=1.2,pt.bg=c("indianred","skyblue3"),bty="n")


#PR(GOT ALOT WORSE)

fit1<-coefs.outL[42:71,1]
fit2<-coefs.outC[42:71,1]
dif<-fit1-fit2



plot(0:29,fit1,ylim=c(0,1),type="l",main=NA,bty="n",
	xlab="Month post outbreak of crisis",ylab="Average Predicted Probability",font=1)
title(main=expression("Pr(Got a lot worse | PartyID)"),font=1)

lines(0:29,fit1,col="indianred",lwd=2)
points(0:29,fit1,bg="indianred",pch=21)
lines(0:29,fit2,col="skyblue3",lwd=2)
points(0:29,fit2,bg="skyblue3",pch=21)


#DIF1
fit1<-coefs.outL[42:71,2]+coefs.outL[42:71,1]
fit2<-coefs.outC[42:71,2]+coefs.outC[42:71,1]
dif<-fit2-fit1



plot(0:29,dif,ylim=c(0,.5),type="l",main=NA,bty="n",
	xlab="Month post outbreak of crisis",ylab="Difference",font=1)
title(main="Partisan Gap",font=1)

lines(0:29,dif,col="black",lwd=2)
points(0:29,dif,bg="black",pch=21)



#DIF2
fit1<-coefs.outL[42:71,1]
fit2<-coefs.outC[42:71,1]
dif<-fit2-fit1



plot(0:29,dif,ylim=c(0,.5),type="l",main=NA,bty="n",
	xlab="Month post outbreak of crisis",ylab="Difference",font=1)
title(main="Partisan Gap",font=1)

lines(0:29,dif,col="black",lwd=2)
points(0:29,dif,bg="black",pch=21)




###################################################################
# APPENDIX J: INCLUDING LIBERAL DEMOCRATS AND NON-IDENTIFIERS	#


#numbers for table A8
#see appendix table A8 for value labels
tab<-table(newdat$party.id)
ptab<-prop.table(tab)





####################
# FIGURE A10	 #

lmfit<-lm(natretro~party.id4*date+party.str+polatt+gender+age+income+mortgage+unemp+education,data=newdat)
summary(lmfit)


#####################
# CONDITIONAL MEANS #

preddat<-data.frame(expand.grid("party.id4"=as.factor(0:3),"date"=as.factor(1:71),
	"party.str"=median(newdat$party.str,na.rm=T),
	"polatt"=median(newdat$polatt,na.rm=T),
	"gender"=median(newdat$gender,na.rm=T),
	"age"=mean(newdat$age,na.rm=T),
	"income"=median(newdat$income,na.rm=T),
	"mortgage"=median(newdat$mortgage,na.rm=T),
	"unemp"=median(newdat$unemp,na.rm=T),
	"education"=median(newdat$education,na.rm=T)
))

preddat.out<-cbind(preddat,fit=predict(lmfit,newdata=preddat,se.fit=T)$fit,
	se=predict(lmfit,newdata=preddat,se.fit=T)$se)
coefs.out<-matrix(nr=71,nc=8)
coefs.out[,1]<-preddat.out$fit[preddat.out$party.id4==0]
#coefs.out[,2]<-preddat.out$se[preddat.out$party.id4==0]
coefs.out[,3]<-preddat.out$fit[preddat.out$party.id4==1]
#coefs.out[,4]<-preddat.out$se[preddat.out$party.id4==1]
coefs.out[,5]<-preddat.out$fit[preddat.out$party.id4==2]
#coefs.out[,6]<-preddat.out$se[preddat.out$party.id4==2]
coefs.out[,7]<-preddat.out$fit[preddat.out$party.id4==3]
#coefs.out[,8]<-preddat.out$se[preddat.out$party.id4==3]

colnames(coefs.out)<-c("Lab","se","Con","se","Libdem","se","Ind","se")


#### PLOT ####

#pdf('retro-all.pdf',width=8,height=5)

par(mfrow=c(1,1),omd=c(0,1,0,1))
plot(coefs.out[,1],ylim=c(0,.8),bty="n",xaxt="n",type="p",lwd=1,xlab="",
	ylab="Predicted Mean",
	col="white")


lines(SMA(coefs.out[,1],n=4),col="red",lwd=3,lty=1)
lines(SMA(coefs.out[,3],n=4),col="blue",lwd=3,lty=1)
lines(SMA(coefs.out[,5],n=4),col="yellow",lwd=3,lty=1)
lines(SMA(coefs.out[,7],n=4),col="darkgrey",lwd=3,lty=1)

points(coefs.out[,1],bg="red",pch=22,cex=.8)
points(coefs.out[,3],bg="blue",pch=23,cex=.8)
points(coefs.out[,5],bg="yellow",pch=24,cex=.8)
points(coefs.out[,7],bg="darkgrey",pch=25,cex=.8)

axis(1,c(1,10,21,32,44,56,68,79,91,97),labels=c("",2005,2006,2007,2008,2009,2010,2011,2012,""),cex.axis=0.8)
text(y=.9,x=40,labels="Northern Rock",cex=1,xpd=TRUE,font=2)
abline(v=40,lwd=2,lty=2)

legend("topright",legend=c("Labour","Conservative","Liberal Democrat","No Identification"),
	lwd=2,
	pt.cex=1,cex=0.8,bty="n",pt.bg=c("red","blue","yellow","darkgrey"),
	col=c("red","blue","yellow","darkgrey"),xpd=TRUE)

points(rep(80,4),c(.79,.745,.7,.655),pch=c(22,23,24,25),bg=c("red","blue","yellow","darkgrey"))


#dev.off()


###########################################
# ATTRIBUTION OF RESPONSIBILITY MEASURE  	#
#	FIGURE A11					#




lmfit<-glm(rpecon2~party.id4*date+party.str+polatt+gender+age+income+mortgage+unemp+education,data=newdat,
	family=binomial(link="probit"))

summary(lmfit)

#####################
# CONDITIONAL MEANS #

preddat.out<-cbind(newdat[rownames(newdat)%in%names(predict(lmfit)),],fit=predict(lmfit,se.fit=T,type="response")$fit,
	se=predict(lmfit,se.fit=T,type="response")$se)


coefs.out<-matrix(NA,nc=4,nr=71) #matrix to hold overall output for plotting

#average predicted probability within groups
for (i in 1:71){
coefs.out[i,1]<-mean(preddat.out$fit[preddat.out$party.id4==0 & preddat.out$date==i])
coefs.out[i,2]<-mean(preddat.out$fit[preddat.out$party.id4==1 & preddat.out$date==i])
coefs.out[i,3]<-mean(preddat.out$fit[preddat.out$party.id4==2 & preddat.out$date==i])
coefs.out[i,4]<-mean(preddat.out$fit[preddat.out$party.id4==3 & preddat.out$date==i])
}

colnames(coefs.out)<-c("Lab","Cons","LibDem","Other")


#### PLOT ####

#pdf('resp-all.pdf',width=8,height=5)

par(mfrow=c(1,1),omd=c(0,1,0,1))
plot(coefs.out[,1],ylim=c(0.2,1),bty="n",xaxt="n",type="p",lwd=1,xlab="",
	ylab="Average Predicted Probability",
	col="white")


lines(SMA(coefs.out[,1],n=4),col="red",lwd=3,lty=1)
lines(SMA(coefs.out[,2],n=4),col="blue",lwd=3,lty=1)
lines(SMA(coefs.out[,3],n=4),col="yellow",lwd=3,lty=1)
lines(SMA(coefs.out[,4],n=4),col="darkgrey",lwd=3,lty=1)

points(coefs.out[,1],bg="red",pch=22,cex=.8)
points(coefs.out[,2],bg="blue",pch=23,cex=.8)
points(coefs.out[,3],bg="yellow",pch=24,cex=.8)
points(coefs.out[,4],bg="darkgrey",pch=25,cex=.8)

axis(1,c(1,10,21,32,44,56,68,79,91,97),labels=c("",2005,2006,2007,2008,2009,2010,2011,2012,""),cex.axis=0.8)
text(y=1.1,x=40,labels="Northern Rock",cex=1,xpd=TRUE,font=2)
#text(y=1.1,x=71,labels="New Government",cex=1,xpd=TRUE,font=2)
abline(v=40,lwd=1,lty=2)
#abline(v=71,lwd=2,lty=2)


legend("topright",legend=c("Labour","Conservative","Liberal Democrat","No Identification"),
	lwd=2,
	pt.cex=1,cex=0.8,bty="n",pt.bg=c("red","blue","yellow","darkgrey"),
	col=c("red","blue","yellow","darkgrey"),xpd=TRUE)

points(rep(80,4),c(.99,.945,.9,.855),pch=c(22,23,24,25),bg=c("red","blue","yellow","darkgrey"))


#dev.off()



###################################################################
# APPENDIX K: CONDITIONING EFFECTS OF THE INFORMATION ENVIRONMENT #


# TABLE A9 #

newdat$readnwspL<-as.factor(newdat$readnwsp)
newdat$readnwspL<-relabel(newdat$readnwspL,
	"1"='Every day',
	"2"='Sometimes',
	"3"='Not at all')

(t<-table(newdat$readnwspL[newdat$date%in%1:71&!is.na(newdat$party.id2)]))
(tp<-round(prop.table(t),digits=2))
sum(t)
require(xtable)
xtable(t)



#	FIGURE A12 	#


#NO EXPOSURE TO NWSP##

lmfit<-lm(natretro~party.id2*date+polatt+gender+age+income+mortgage+unemp+education
	,data=newdat[newdat$readnwsp==3,])
summary(lmfit)




preddat<-data.frame(expand.grid("party.id2"=c(0,1),"date"=as.factor(1:71),
	"party.str"=median(newdat$party.str,na.rm=T),
	"polatt"=median(newdat$polatt,na.rm=T),
	"gender"=median(newdat$gender,na.rm=T),
	"age"=mean(newdat$age,na.rm=T),
	"income"=median(newdat$income,na.rm=T),
	"mortgage"=median(newdat$mortgage,na.rm=T),
	"unemp"=median(newdat$unemp,na.rm=T),
	"education"=median(newdat$education,na.rm=T)
))

vmat<-vcov(lmfit)


preddat.out<-cbind(preddat,predict.new(lmfit,vmat,preddat))


coefs.out<-matrix(NA,nc=6,nr=71) #matrix to hold overall output for plotting
coefs.out[,1]<-preddat.out$fit[preddat.out$party.id2==0]
coefs.out[,2]<-preddat.out$se.fit[preddat.out$party.id2==0]
coefs.out[,3]<-preddat.out$fit[preddat.out$party.id2==1]
coefs.out[,4]<-preddat.out$se.fit[preddat.out$party.id2==1]

# PARTISAN GAP #

est<-matrix(c(coef(lmfit)["party.id2"],
	rep(coef(lmfit)["party.id2"],70)+
	coef(lmfit)[grep("^party.id2:date",names(coef(lmfit)))]
		))


# SEs for marginal effect #


b1<-rep(vmat["party.id2","party.id2"],70)
b3<-diag(vmat)[grep("^party.id2:date",rownames(vmat))]
b1b3<-vmat["party.id2",grep("^party.id2:date",rownames(vmat))] #extract covar

se<-cbind(b1,b3,b1b3);rownames(se)<-rep("",70) #arrange in matrix form

se.out<-apply(se,1,function(x) sqrt(x[1]+x[2]+2*x[3]))
se.out<-c(b1[1],se.out)

# Final output #

est.out<-cbind(est,se.out)
coefs.out[,5:6]<-est.out;coefs.out<-round(coefs.out,digits=5)
coefs.out[,5]<- -coefs.out[,5]
colnames(coefs.out)<-c("LAB","SE","CON","SE","GAP","SE")


coefsWEAK<-coefs.out 	#storing





######MIDDLING IDENTIFIERS########



lmfit<-lm(natretro~party.id2*date+polatt+gender+age+income+mortgage+unemp+education
	,data=newdat[newdat$readnwsp==2,])
summary(lmfit)




#####################
# CONDITIONAL MEANS #

preddat<-data.frame(expand.grid("party.id2"=c(0,1),"date"=as.factor(1:71),
	"party.str"=median(newdat$party.str,na.rm=T),
	"polatt"=median(newdat$polatt,na.rm=T),
	"gender"=median(newdat$gender,na.rm=T),
	"age"=mean(newdat$age,na.rm=T),
	"income"=median(newdat$income,na.rm=T),
	"mortgage"=median(newdat$mortgage,na.rm=T),
	"unemp"=median(newdat$unemp,na.rm=T),
	"education"=median(newdat$education,na.rm=T)
))

vmat<-vcov(lmfit)


preddat.out<-cbind(preddat,predict.new(lmfit,vmat,preddat))


coefs.out<-matrix(NA,nc=6,nr=71) #matrix to hold overall output for plotting
coefs.out[,1]<-preddat.out$fit[preddat.out$party.id2==0]
coefs.out[,2]<-preddat.out$se.fit[preddat.out$party.id2==0]
coefs.out[,3]<-preddat.out$fit[preddat.out$party.id2==1]
coefs.out[,4]<-preddat.out$se.fit[preddat.out$party.id2==1]

################
# PARTISAN GAP #

est<-matrix(c(coef(lmfit)["party.id2"],
	rep(coef(lmfit)["party.id2"],70)+
	coef(lmfit)[grep("^party.id2:date",names(coef(lmfit)))]
		))


###########################
# SEs for marginal effect #


b1<-rep(vmat["party.id2","party.id2"],70)
b3<-diag(vmat)[grep("^party.id2:date",rownames(vmat))]
b1b3<-vmat["party.id2",grep("^party.id2:date",rownames(vmat))] #extract covar

se<-cbind(b1,b3,b1b3);rownames(se)<-rep("",70) #arrange in matrix form

se.out<-apply(se,1,function(x) sqrt(x[1]+x[2]+2*x[3]))
se.out<-c(b1[1],se.out)

################
# Final output #

est.out<-cbind(est,se.out)
coefs.out[,5:6]<-est.out;coefs.out<-round(coefs.out,digits=5)
coefs.out[,5]<- -coefs.out[,5]
colnames(coefs.out)<-c("LAB","SE","CON","SE","GAP","SE")

coefsMID<-coefs.out





###############################
# 	HIGH EXPOSURE			#


lmfit<-lm(natretro~party.id2*date+polatt+gender+age+income+mortgage+unemp+education
	,data=newdat[newdat$readnwsp==1,])
#summary(lmfit)




#####################
# CONDITIONAL MEANS #

preddat<-data.frame(expand.grid("party.id2"=c(0,1),"date"=as.factor(1:71),
	"party.str"=median(newdat$party.str,na.rm=T),
	"polatt"=median(newdat$polatt,na.rm=T),
	"gender"=median(newdat$gender,na.rm=T),
	"age"=mean(newdat$age,na.rm=T),
	"income"=median(newdat$income,na.rm=T),
	"mortgage"=median(newdat$mortgage,na.rm=T),
	"unemp"=median(newdat$unemp,na.rm=T),
	"education"=median(newdat$education,na.rm=T)
))

vmat<-vcov(lmfit)


preddat.out<-cbind(preddat,predict.new(lmfit,vmat,preddat))


coefs.out<-matrix(NA,nc=6,nr=71) #matrix to hold overall output for plotting
coefs.out[,1]<-preddat.out$fit[preddat.out$party.id2==0]
coefs.out[,2]<-preddat.out$se.fit[preddat.out$party.id2==0]
coefs.out[,3]<-preddat.out$fit[preddat.out$party.id2==1]
coefs.out[,4]<-preddat.out$se.fit[preddat.out$party.id2==1]

################
# PARTISAN GAP #

est<-matrix(c(coef(lmfit)["party.id2"],
	rep(coef(lmfit)["party.id2"],70)+
	coef(lmfit)[grep("^party.id2:date",names(coef(lmfit)))]
		))


###########################
# SEs for marginal effect #


b1<-rep(vmat["party.id2","party.id2"],70)
b3<-diag(vmat)[grep("^party.id2:date",rownames(vmat))]
b1b3<-vmat["party.id2",grep("^party.id2:date",rownames(vmat))] #extract covar

se<-cbind(b1,b3,b1b3);rownames(se)<-rep("",70) #arrange in matrix form

se.out<-apply(se,1,function(x) sqrt(x[1]+x[2]+2*x[3]))
se.out<-c(b1[1],se.out)

# Final output #

est.out<-cbind(est,se.out)
coefs.out[,5:6]<-est.out;coefs.out<-round(coefs.out,digits=5)
coefs.out[,5]<- -coefs.out[,5]
colnames(coefs.out)<-c("LAB","SE","CON","SE","GAP","SE")
coefsSTRONG<-coefs.out


# PLOTTING	#


#pdf('natretro-nwspexp-composite.pdf',width=10,height=6.5)


col1<-rgb(.2,.2,.2,alpha=.8)
col2<-rgb(.4,.4,.4,alpha=.8)
d1<-.4
d2<-.3
d3<-.6
d4<-.7

td<-1.75
subd<-0.25
par(mfrow=c(1,1))
par(omd=c(0,d1,.375,1))
coefs.out<-coefsWEAK
plot(coefs.out[,1],ylim=c(0,1.15),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="Predicted Mean",
	col=col1)
axis(1,labels=NA)
mtext(3,text='A',line=td,font=2,cex=1.5)
mtext(3,text='No exposure to newspapers',line=subd,font=1,cex=1)
segments(1:71,coefs.out[,1],1:71,coefs.out[,1]+1.97*coefs.out[,2],col=col1);segments(1:71,coefs.out[,1],1:71,coefs.out[,1]-1.97*coefs.out[,2],col=col1)
segments(1:71,coefs.out[,3],1:71,coefs.out[,3]+1.97*coefs.out[,4],col=col1);segments(1:71,coefs.out[,3],1:71,coefs.out[,3]-1.97*coefs.out[,4],col=col1)
points(coefs.out[,1],pch=21,bg=col1,col=col1)
points(coefs.out[,3],pch=21,col=col1)
lines(SMA(coefs.out[,1],n=2))
lines(SMA(coefs.out[,3],n=2),lty=1)
abline(h=1.05)
mtext(3,text="Retrospective Economic Perceptions",line=-1,font=1,cex=.8)
legend(-3,1.05,pch=c(21,1),pt.bg=c(col1,NA),col=col1,legend=c("Labour Identifier","Conservative Identifier"),cex=.9,bty="n")
par(omd=c(0,d1,.375,.95),new=TRUE)
abline(v=40,lty=2)
par(omd=c(0,d1,.375,1),new=TRUE)
par(omd=c(0,d1,0,.625),new=TRUE)
plot(coefs.out[,5],ylim=c(-.1,.6),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="Marginal Effect of Party Id.",yaxt='n',
	col="white")
axis(2,at=seq(-.4,.4,.2),las=2)
abline(h=0,lty=2)
abline(h=.55)
mtext(3,text="Partisan Differences",line=-1,font=1,cex=.8)
par(omd=c(0,d1,0,.575),new=TRUE)
abline(v=40,lty=2)
par(omd=c(0,d1,0,.625),new=TRUE)
segments(1:71,coefs.out[,5],1:71,coefs.out[,5]+1.97*coefs.out[,6],col=col2);segments(1:71,coefs.out[,5],1:71,coefs.out[,5]-1.97*coefs.out[,6],col=col2)
lines(SMA(coefs.out[,5],n=2))
points(coefs.out[,5],pch=21,bg=col2,col=col2)
axis(1,c(10,21,32,44,56,68,79,91),labels=c(2005,"",2007,"",2009,"",2011,""),cex.axis=1)

##ADDING MIDDLING CATEGORY###
par(omd=c(d2,d4,.375,1),new=T)
coefs.out<-coefsMID
plot(coefs.out[,1],ylim=c(0,1.15),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="",yaxt="n",
	col=col1)
axis(1,labels=NA)
axis(2,labels=NA)
mtext(3,text='B',line=td,font=2,cex=1.5)
mtext(3,text='Middling exposure to newspapers',line=subd,font=1,cex=1)
segments(1:71,coefs.out[,1],1:71,coefs.out[,1]+1.97*coefs.out[,2],col=col1);segments(1:71,coefs.out[,1],1:71,coefs.out[,1]-1.97*coefs.out[,2],col=col1)
segments(1:71,coefs.out[,3],1:71,coefs.out[,3]+1.97*coefs.out[,4],col=col1);segments(1:71,coefs.out[,3],1:71,coefs.out[,3]-1.97*coefs.out[,4],col=col1)
points(coefs.out[,1],pch=21,bg=col1,col=col1)
points(coefs.out[,3],pch=21,col=col1)
lines(SMA(coefs.out[,1],n=2))
lines(SMA(coefs.out[,3],n=2),lty=1)
abline(h=1.05)
mtext(3,text="Retrospective Economic Perceptions",line=-1,font=1,cex=.8)
text(0,.9,labels="Northern Rock",font=1,cex=.8,pos=4)
arrows(30,.875,40,.8,length=.075)
par(omd=c(d2,d4,.375,.95),new=TRUE)
abline(v=40,lty=2)
par(omd=c(d2,d4,.375,1),new=TRUE)
par(omd=c(d2,d4,0,.625),new=TRUE)
plot(coefs.out[,5],ylim=c(-.1,.6),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="",yaxt='n',
	col="white")
axis(2,at=seq(-.4,.4,.2),labels=NA,las=2)
abline(h=0,lty=2)
abline(h=.55)
mtext(3,text="Partisan Differences",line=-1,font=1,cex=.8)
par(omd=c(d2,d4,0,.575),new=TRUE)
abline(v=40,lty=2)
par(omd=c(d2,d4,0,.625),new=TRUE)
segments(1:71,coefs.out[,5],1:71,coefs.out[,5]+1.97*coefs.out[,6],col=col2);segments(1:71,coefs.out[,5],1:71,coefs.out[,5]-1.97*coefs.out[,6],col=col2)
lines(SMA(coefs.out[,5],n=2))
points(coefs.out[,5],pch=21,bg=col2,col=col2)
axis(1,c(10,21,32,44,56,68,79,91),labels=c(2005,"",2007,"",2009,"",2011,""),cex.axis=1)

##ADDING STRONG CATEGORY###
par(omd=c(d3,1,.375,1),new=T)
coefs.out<-coefsSTRONG
plot(coefs.out[,1],ylim=c(0,1.15),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="",yaxt="n",
	col=col1)
axis(1,labels=NA)
axis(2,labels=NA)
mtext(3,text='C',line=td,font=2,cex=1.5)
mtext(3,text='High exposure to newspapers',line=subd,font=1,cex=1)
segments(1:71,coefs.out[,1],1:71,coefs.out[,1]+1.97*coefs.out[,2],col=col1);segments(1:71,coefs.out[,1],1:71,coefs.out[,1]-1.97*coefs.out[,2],col=col1)
segments(1:71,coefs.out[,3],1:71,coefs.out[,3]+1.97*coefs.out[,4],col=col1);segments(1:71,coefs.out[,3],1:71,coefs.out[,3]-1.97*coefs.out[,4],col=col1)
points(coefs.out[,1],pch=21,bg=col1,col=col1)
points(coefs.out[,3],pch=21,col=col1)
lines(SMA(coefs.out[,1],n=2))
lines(SMA(coefs.out[,3],n=2),lty=1)
abline(h=1.05)
mtext(3,text="Retrospective Economic Perceptions",line=-1,font=1,cex=.8)
par(omd=c(d3,1,.375,.95),new=TRUE)
abline(v=40,lty=2)
par(omd=c(d3,1,.375,1),new=TRUE)
par(omd=c(d3,1,0,.625),new=TRUE)
plot(coefs.out[,5],ylim=c(-.1,.6),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="",yaxt='n',
	col="white")
axis(2,at=seq(-.4,.4,.2),labels=NA,las=2)
abline(h=0,lty=2)
abline(h=.55)
mtext(3,text="Partisan Differences",line=-1,font=1,cex=.8)
par(omd=c(d3,1,0,.575),new=TRUE)
abline(v=40,lty=2)
par(omd=c(d3,1,0,.625),new=TRUE)
segments(1:71,coefs.out[,5],1:71,coefs.out[,5]+1.97*coefs.out[,6],col=col2);segments(1:71,coefs.out[,5],1:71,coefs.out[,5]-1.97*coefs.out[,6],col=col2)
lines(SMA(coefs.out[,5],n=2))
points(coefs.out[,5],pch=21,bg=col2,col=col2)
axis(1,c(10,21,32,44,56,68,79,91),labels=c(2005,"",2007,"",2009,"",2011,""),cex.axis=1)

#dev.off()





###################################
#		ATTRIBUTIONS 	####

boot=100 #runs slow with 1000 iterations

#collapse date into 6m bins
newdat$datehy<-recode(as.numeric(newdat$date),
	'1:6=0;
	7:12=1;
	13:17=2; 
	18:23=3;
	24:28=4; 
	29:34=5;
	35:40=6;
	41:46=7;
	47:52=8;
	53:58=9;
	59:64=10;
	65:71=11;else=NA')
table(newdat$date,newdat$datehy)#check joint distribution with date
newdat$datehy<-factor(newdat$datehy)



#  NO EXP TO NWSP  #

lmfit<-glm(rpecon2~party.id2*datehy+polatt+coneu+gender+age+income+mortgage+unemp+education,
	data=newdat[newdat$readnwsp==3,],
	family=binomial(link="probit"))

mtable(lmfit)

coefsrp<-matrix(NA,nc=6,nr=12) #matrix to hold overall output for plotting

varnam<-c("coneu","party.str","readnwsp","mortgage","education","income","gender","age","polatt","unemp","party.id2","datehy")

# APP and AME  #

	x<-lmfit 	#Set model
	vmat<-vcov(x)
	samp<-mvrnorm(n=boot,mu=coef(x),Sigma=vmat) #simulate estimation uncertainty
dim(vmat)
	#set up datehy frame with observered values
	datx1<-newdat[complete.cases(newdat[,varnam]),varnam];datx1$party.id2<-1 #data where all x's are set to 1
	datx0<-newdat[complete.cases(newdat[,varnam]),varnam];datx0$party.id2<-0 #set to 0

	ame.samp<-rep(NA,boot)
	app1<-rep(NA,boot)
	app0<-rep(NA,boot)
	ame<-rep(NA,boot)

	for (i in 1:12){
	for (k in 1:boot){
		x$coefficients<-samp[k,]
		app0[k]<-mean(predict(x,newdat=datx0[datx0$datehy==i-1,],type="response"))
		app1[k]<-mean(predict(x,newdat=datx1[datx1$datehy==i-1,],type="response"))
		ame[k]<-app0[k]-app1[k]
	}
		coefsrp[i,]<-c(mean(app0),sd(app0),mean(app1),sd(app1),mean(ame),sd(ame)) #storing time varying APPs and AMEs
	}

coefsrp
coefsrpWEAK<-coefsrp


#  MIDDLING EXP  #

lmfit<-glm(rpecon2~party.id2*datehy+polatt+coneu+gender+age+income+mortgage+unemp+education,
	data=newdat[newdat$readnwsp==2,],
	family=binomial(link="probit"))

mtable(lmfit)

coefsrp<-matrix(NA,nc=6,nr=12) #matrix to hold overall output for plotting

varnam<-c("coneu","party.str","readnwsp","mortgage","education","income","gender","age","polatt","unemp","party.id2","datehy")

# APP and AME  #

	x<-lmfit 	#Set model
	vmat<-vcov(x)
	samp<-mvrnorm(n=boot,mu=coef(x),Sigma=vmat) #simulate estimation uncertainty
dim(vmat)
	#set up datehy frame with observered values
	datx1<-newdat[complete.cases(newdat[,varnam]),varnam];datx1$party.id2<-1 #data where all x's are set to 1
	datx0<-newdat[complete.cases(newdat[,varnam]),varnam];datx0$party.id2<-0 #set to 0

	ame.samp<-rep(NA,boot)
	app1<-rep(NA,boot)
	app0<-rep(NA,boot)
	ame<-rep(NA,boot)

	for (i in 1:12){
	for (k in 1:boot){
		x$coefficients<-samp[k,]
		app0[k]<-mean(predict(x,newdat=datx0[datx0$datehy==i-1,],type="response"))
		app1[k]<-mean(predict(x,newdat=datx1[datx1$datehy==i-1,],type="response"))
		ame[k]<-app0[k]-app1[k]
	}
		coefsrp[i,]<-c(mean(app0),sd(app0),mean(app1),sd(app1),mean(ame),sd(ame)) #storing time varying APPs and AMEs
	}

coefsrp
coefsrpMID<-coefsrp


#  HIGH EXP  #

lmfit<-glm(rpecon2~party.id2*datehy+polatt+coneu+gender+age+income+mortgage+unemp+education,
	data=newdat[newdat$readnwsp==1,],
	family=binomial(link="probit"))

mtable(lmfit)
coefsrp<-matrix(NA,nc=6,nr=12) #matrix to hold overall output for plotting

varnam<-c("coneu","party.str","readnwsp","mortgage","education","income","gender","age","polatt","unemp","party.id2","datehy")

# APP and AME  #

	x<-lmfit 	#Set model
	vmat<-vcov(x)
	samp<-mvrnorm(n=boot,mu=coef(x),Sigma=vmat) #simulate estimation uncertainty
dim(vmat)
	#set up datehy frame with observered values
	datx1<-newdat[complete.cases(newdat[,varnam]),varnam];datx1$party.id2<-1 #data where all x's are set to 1
	datx0<-newdat[complete.cases(newdat[,varnam]),varnam];datx0$party.id2<-0 #set to 0

	ame.samp<-rep(NA,boot)
	app1<-rep(NA,boot)
	app0<-rep(NA,boot)
	ame<-rep(NA,boot)

	for (i in 1:12){
	for (k in 1:boot){
		x$coefficients<-samp[k,]
		app0[k]<-mean(predict(x,newdat=datx0[datx0$datehy==i-1,],type="response"))
		app1[k]<-mean(predict(x,newdat=datx1[datx1$datehy==i-1,],type="response"))
		ame[k]<-app0[k]-app1[k]
	}
		coefsrp[i,]<-c(mean(app0),sd(app0),mean(app1),sd(app1),mean(ame),sd(ame)) #storing time varying APPs and AMEs
	}

coefsrp
coefsrpSTRONG<-coefsrp



#plotting#

coefsrpWEAK[,5]<- -coefsrpWEAK[,5]
coefsrpMID[,5]<- -coefsrpMID[,5]
coefsrpSTRONG[,5]<- -coefsrpSTRONG[,5] #new sign for difference


#pdf('rpecon-nwspexp-composite.pdf',width=10,height=6.5)


col1<-rgb(.2,.2,.2,alpha=.8)
col2<-rgb(.4,.4,.4,alpha=.8)

d1<-.4
d2<-.3
d3<-.6
d4<-.7

td<-1.75
subd<-0.25

par(mfrow=c(1,1))
par(omd=c(0,d1,.375,1))

coefsrp<-coefsrpWEAK

plot(coefsrp[,1],ylim=c(0,1.15),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="Avg. Predicted Probability",
	col=col1)

mtext(3,text='A',line=td,font=2,cex=1.5)
mtext(3,text='No exposure to newspapers',line=subd,font=1,cex=1)

segments(1:11,coefsrp[,1],1:11,coefsrp[,1]+1.97*coefsrp[,2],col=col1);segments(1:11,coefsrp[,1],1:11,coefsrp[,1]-1.97*coefsrp[,2],col=col1)
segments(1:11,coefsrp[,3],1:11,coefsrp[,3]+1.97*coefsrp[,4],col=col1);segments(1:11,coefsrp[,3],1:11,coefsrp[,3]-1.97*coefsrp[,4],col=col1)

points(coefsrp[,1],pch=21,bg=col1,col=col1)
points(coefsrp[,3],pch=21,col=col1)

lines(coefsrp[,1])
lines(coefsrp[,3],lty=1)

adjtick<-c(10,21,32,44,56,68,79,91)/6
axis(1,adjtick,labels=NA)

abline(h=1.05)
mtext(3,text="Attribution of Responsibility",line=-1,font=1,cex=.8)

legend("bottomleft",pch=c(21,1),pt.bg=c(col1,NA),col=col1,legend=c("Labour Identifier","Conservative Identifier"),cex=.9,bty="n")

par(omd=c(0,d1,.375,.95),new=TRUE)
abline(v=40/6,lty=2)
par(omd=c(0,d1,.375,1),new=TRUE)

text(6,.975,labels="Northern Rock",font=1,cex=.8,pos=2,xpd=T)
arrows(4,.925,6.5,.8,length=.075)



par(omd=c(0,d1,0,.625),new=TRUE)

plot(coefsrp[,5],ylim=c(-.5,.65),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="Avg. Marginal Effect",yaxt="n",
	col=col2)
axis(2,at=seq(-.4,.4,.2),las=2)

abline(h=.55)
mtext(3,text="Partisan Differences",line=-1,font=1,cex=.8)

segments(1:71,coefsrp[,5],1:71,coefsrp[,5]+1.97*coefsrp[,6],col=col2);segments(1:71,coefsrp[,5],1:71,coefsrp[,5]-1.97*coefsrp[,6],col=col2)
lines(coefsrp[,5])
points(coefsrp[,5],pch=21,bg=col2,col=col2)
abline(h=0,lty=2)

par(omd=c(0,d1,0,.575),new=TRUE)
abline(v=40/6,lty=2)
par(omd=c(0,d1,0,.625),new=TRUE)


axis(1,adjtick,labels=c(2005,"",2007,"",2009,"",2011,""),cex.axis=1)


########MIDDLING IDers#######
coefsrp<-coefsrpMID

par(omd=c(d2,d4,.375,1),new=TRUE)
plot(coefsrp[,1],ylim=c(0,1.15),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="",yaxt="n",
	col=col1,add=TRUE)

mtext(3,text='B',line=td,font=2,cex=1.5)
mtext(3,text='Middling exposure to newspapers',line=subd,font=1,cex=1)

segments(1:11,coefsrp[,1],1:11,coefsrp[,1]+1.97*coefsrp[,2],col=col1);segments(1:11,coefsrp[,1],1:11,coefsrp[,1]-1.97*coefsrp[,2],col=col1)
segments(1:11,coefsrp[,3],1:11,coefsrp[,3]+1.97*coefsrp[,4],col=col1);segments(1:11,coefsrp[,3],1:11,coefsrp[,3]-1.97*coefsrp[,4],col=col1)

points(coefsrp[,1],pch=21,bg=col1,col=col1)
points(coefsrp[,3],pch=21,col=col1)

lines(coefsrp[,1])
lines(coefsrp[,3],lty=1)

adjtick<-c(10,21,32,44,56,68,79,91)/6
axis(1,adjtick,labels=NA)
axis(2,labels=NA)

abline(h=1.05)
mtext(3,text="Attribution of Responsibility",line=-1,font=1,cex=.8)


par(omd=c(d2,d4,.375,.95),new=TRUE)
abline(v=40/6,lty=2)
par(omd=c(d2,d4,.375,1),new=TRUE)


par(omd=c(d2,d4,0,.625),new=TRUE)

plot(coefsrp[,5],ylim=c(-.5,.65),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="",yaxt="n",
	col=col2)
axis(2,at=seq(-.4,.4,.2),labels=NA,las=2)

abline(h=.55)
mtext(3,text="Partisan Differences",line=-1,font=1,cex=.8)

segments(1:71,coefsrp[,5],1:71,coefsrp[,5]+1.97*coefsrp[,6],col=col2);segments(1:71,coefsrp[,5],1:71,coefsrp[,5]-1.97*coefsrp[,6],col=col2)
lines(coefsrp[,5])
points(coefsrp[,5],pch=21,bg=col2,col=col2)
abline(h=0,lty=2)

par(omd=c(d2,d4,0,.575),new=TRUE)
abline(v=40/6,lty=2)
par(omd=c(d2,d4,0,.625),new=TRUE)


axis(1,adjtick,labels=c(2005,"",2007,"",2009,"",2011,""),cex.axis=1)



####### 	ADDING STRONG ID GRAPH #############

coefsrp<-coefsrpSTRONG

par(omd=c(d3,1,.375,1),new=TRUE)
plot(coefsrp[,1],ylim=c(0,1.15),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="",yaxt="n",
	col=col1,add=TRUE)

mtext(3,text='C',line=td,font=2,cex=1.5)
mtext(3,text='High exposure to newspapers',line=subd,font=1,cex=1)

segments(1:11,coefsrp[,1],1:11,coefsrp[,1]+1.97*coefsrp[,2],col=col1);segments(1:11,coefsrp[,1],1:11,coefsrp[,1]-1.97*coefsrp[,2],col=col1)
segments(1:11,coefsrp[,3],1:11,coefsrp[,3]+1.97*coefsrp[,4],col=col1);segments(1:11,coefsrp[,3],1:11,coefsrp[,3]-1.97*coefsrp[,4],col=col1)

points(coefsrp[,1],pch=21,bg=col1,col=col1)
points(coefsrp[,3],pch=21,col=col1)

lines(coefsrp[,1])
lines(coefsrp[,3],lty=1)

adjtick<-c(10,21,32,44,56,68,79,91)/6
axis(1,adjtick,labels=NA)
axis(2,labels=NA)

abline(h=1.05)
mtext(3,text="Attribution of Responsibility",line=-1,font=1,cex=.8)


par(omd=c(d3,1,.375,.95),new=TRUE)
abline(v=40/6,lty=2)
par(omd=c(d3,1,.375,1),new=TRUE)


par(omd=c(d3,1,0,.625),new=TRUE)

plot(coefsrp[,5],ylim=c(-.5,.65),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="",yaxt="n",
	col=col2)
axis(2,at=seq(-.4,.4,.2),labels=NA,las=2)

abline(h=.55)
mtext(3,text="Partisan Differences",line=-1,font=1,cex=.8)

segments(1:71,coefsrp[,5],1:71,coefsrp[,5]+1.97*coefsrp[,6],col=col2);segments(1:71,coefsrp[,5],1:71,coefsrp[,5]-1.97*coefsrp[,6],col=col2)
lines(coefsrp[,5])
points(coefsrp[,5],pch=21,bg=col2,col=col2)
abline(h=0,lty=2)

par(omd=c(d3,1,0,.575),new=TRUE)
abline(v=40/6,lty=2)
par(omd=c(d3,1,0,.625),new=TRUE)


axis(1,adjtick,labels=c(2005,"",2007,"",2009,"",2011,""),cex.axis=1)


#dev.off()


###################################################
# APPENDIX L: Full model output for main analysis #
# see main script for the model output.		  # 



